function [loglike,TCPmodel,EQD2,Effect,ps] = loglikelihood_logistic(par_long,data_TCPexp,inds,dep_a,dep_b)%[loglike,TCPmodel,EQD2] = loglikelihood_logistic(par_long,data_TCPexp)% n,d,pacs,TCPexp,N,param)
%This function calculates the fitting loglikelihood of a logistic TCP model using some input parameters (par_long and N) and a 
%SF expresssion with functional dependencies dep_a, and dep_b for every cohort of a clinical treatment database (data_TCPexp).

%Input:

%        1. par_long: TCP model initial parameter values par_long=[alpha/beta a b D50 gamma lambda_prolif t_kickoff TCPmax]
%        2. data_TCPexp: database array with information of clinical RT treatment for several patient cohorts 
%        3. inds: Cohorts in the database with more than 170 pts.
%        4. dep_a: functional dependence with dose of the alpha term in the expression of the surviving fraction
%        5. dep_b: functional dependence with dose of the beta term in the expression of the surviving fraction



%Output: 
%        1. loglike: loglikelihood of the best fit
%        2. TCPmodel: TCP model values associated to each cohort of the database calcualted using dep_a and dep_b model functional dependencies and the input parameter values         
%        3. EQD2: EQD associated to each cohort of the database calculated using dep_a and dep_b model functional dependencies and the input parameter values
%        4. Effect_opt: Effect,log(SF), associated to each cohort of the database, calculated using dep_a and dep_b model functional dependencies and the input parameter values         
%        5. ps: fitting probability (use of binomial distribution) 


%Model parameters:
alphabeta=par_long(1);
a=par_long(2);
b=par_long(3);
D50=par_long(4);
gamma_model=par_long(5);
lambda=par_long(6);
tk=par_long(7);
TCPmax=par_long(8);
%Assumed alpha radiosensitivity value
alpha=0.3; 
%Treatment information from every cohort in the clinical database: number of treatment schedules, number of fractions, dose per fraction, number of patients/metastases involved, treatment time and tumor control probability. 
schedules=data_TCPexp(1,7); %number of treatment schedules
n=data_TCPexp(:,7+1:2:7+(2*schedules-1)); %number of fractions
d=data_TCPexp(:,7+2:2:7+2*schedules); %average dose per fratcion
pacs=data_TCPexp(:,1); %number of treated patients for NSCLC/ treated metastates for brain
pTCPexp=data_TCPexp(:,3); %Agresti-Coul probability associated to experimental Kaplan TCP value
t_trat=data_TCPexp(:,6); %Treatment time for calculation of accelerated repopulation effect
num_datapoints=length(data_TCPexp(:,7));
vector_ones=ones(max(data_TCPexp(:,1)),1);

%Calculation of treatment effects and EQD2 values
Effect=max(sum(n.*d.*alpha.*(1+a.*d.^dep_a+d.*(1+b.*d.^dep_b)/alphabeta),2)-lambda*max(0,t_trat-tk),0);
    if dep_a==-1 %LQL model
         Effect=max(sum(n.*alpha.*(d+2.*(a.*d+exp(-a.*d)-1)./(alphabeta.*a^2)),2)-lambda*max(0,t_trat-tk),0);  
    end  
[EQD2]=EQD2_calculator(Effect,[alphabeta a b 0 0 lambda tk alpha],dep_a,dep_b);
    if any(isnan(EQD2))==false %if no NaN, calculate TCP, prob, loglike..
         TCPmodel=TCPmax./(1+(D50./EQD2').^(4.*gamma_model));
         %Calculation of the probability that the experimental TCP value yielded the number of controls estimated from the radiobiological
         %model TCPmodel. Use of binomial probability, where factor is all x (tumor control cases)-combinations of a set Ns patients or metastases.
         factor=gamma(pacs+1)./(gamma(pacs.*TCPmodel+1).*(gamma(pacs-pacs.*TCPmodel+1)));
         %if controls=0, factor=1
         inds_with_0controls=inds(round(TCPmodel(inds).*pacs(inds))==0);
         factor(inds_with_0controls)=1;
         %delete inds for those patients
         idx=ismember(inds,inds_with_0controls);
         inds(idx)=[];
         %if control=1, factor=pacs
         inds_with_1control=inds(round(TCPmodel(inds).*pacs(inds))==1);
         factor(inds_with_1control)=pacs(inds_with_1control);
         idx=ismember(inds,inds_with_1control);
         inds(idx)=[];
         %binomial fit probability calculation  
         ps=factor.*pTCPexp.^(pacs.*TCPmodel).*(1-pTCPexp).^(pacs.*(1-TCPmodel));
         %To avoid divergence (Inf) when patient/metastases number (Ns) is very large, round x (with small effect) and calculate 
         %factor=n!/(x!(n-x)!)=prod([max(x,n-x)+1:1:n))/factorial(min(x,n-x)) with a=max(x,n-x)and b=round(min(x,n-x)); 
             
         for ij=1:length(inds) 
                   xr=pacs(inds(ij)).*TCPmodel(inds(ij));
                   x=fix(xr);  
                   %If any of the values goes to infinite, use binopdf
                   % first interpolation point (x-1), avoiding negative values
                   x0=max(0,x-1);
                   pacs_i=pacs(inds(ij));
                   pTCPexp_i=pTCPexp(inds(ij));
                   n_mx=pacs_i-x0;
                   a_f=max(x0,n_mx)+1;
                   b_f=min(x0,n_mx);
                   vector_j=a_f:1:pacs_i;
                   vector_num=vector_ones;
                   vector_num(1:length(vector_j))=vector_j;
                   vector_den=vector_ones;
                   vector_den(1:b_f)=1:1:b_f;
                   factor_1=prod(vector_num./vector_den);
                   p1=factor_1.*pTCPexp_i.^(x0).*(1-pTCPexp_i).^(n_mx);

                   if p1<Inf
                   %second interpolation point (x), central point (c)
                      n_mx=pacs_i-x;
                      a_f=max(x,n_mx)+1;
                      b_f=min(x,n_mx);
                      vector_j=a_f:1:pacs_i;
                      vector_num=vector_ones;
                      vector_num(1:length(vector_j))=vector_j;
                      vector_den=vector_ones;
                      vector_den(1:b_f)=1:1:b_f;
                      factor_c=prod(vector_num./vector_den); 
                      pc=factor_c.*pTCPexp_i.^(x).*(1-pTCPexp_i).^(n_mx);

                      if pc<Inf
                         % third interpolation point (x+1), avoiding x+1>pacs 
                         x1=min(pacs_i,x+1);
                         n_mx=pacs_i-x1;
                         a_f=max(x1,n_mx)+1;
                         b_f=min(x1,n_mx);
                         vector_j=a_f:1:pacs_i;
                         vector_num=vector_ones;
                         vector_num(1:length(vector_j))=vector_j;
                         vector_den=vector_ones;
                         vector_den(1:b_f)=1:1:b_f;
                         factor_2=prod(vector_num./vector_den);
                         p2=factor_2.*pTCPexp_i.^(x1).*(1-pTCPexp_i).^(n_mx);

                         if p2< Inf
                             p_v=[p1 pc p2];
                             %Avoid repeated values
                             [x_unique,i_unique]=unique([x0 x x1]);
                             p_unique=p_v(i_unique);
                             ps(inds(ij))=spline(x_unique,p_unique,xr);
                              if ps(inds(ij))< 0
                                 ps(inds(ij))=spline([x x1],[pc p2],xr);
                              end
                         else
                          ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                         end
                      else 
                       ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                      end
                   else 
                     ps(inds(ij))=calc_binopdf(xr,pacs_i,pTCPexp_i);
                   end  
             end
         loglike=-log(prod(ps));
    else %If there is no EQD2, discard parameters (loglike set to large value, and variables = NaN or 0)
         TCPmodel=NaN(num_datapoints,1);
         EQD2=NaN(num_datapoints,1);
         loglike=1e100;
         ps=zeros(num_datapoints,1);
         Effect=NaN(num_datapoints,1);
    end
end

function [p]=calc_binopdf(xr,n,p)
x=fix(xr);
x0=max(0,x-1);
x1=min(n,x+1);
p1=binopdf(x0,n,p);
pc=binopdf(x,n,p);
p2=binopdf(x1,n,p);
p_v=[p1 pc p2];
%Avoid repeated values
[x_unique,i_unique]=unique([x0 x x1]);
p_unique=p_v(i_unique);
p=spline(x_unique,p_unique,xr);
%Avoid interpolation of a negative value
if p< 0
    p=spline([x x1],[pc p2],xr);
end
end
    